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Abstract 

Molecular dynamics simulations are performed to investigate the plastic response of a model 
glass to a local shear transformation in a quiescent system. The deformation of the material is 
induced by a spherical inclusion that is gradually strained into an ellipsoid of the same volume and 
then reverted back into the sphere. We show that the number of cage-breaking events increases 
with increasing strain amplitude of the shear transformation. The results of numerical simulations 
indicate that the density of cage jumps is larger in the cases of weak damping or slow shear 
transformation. Remarkably, we also found that, for a given strain amplitude, the peak value of 
the density profiles is a function of the ratio of the damping coefficient and the time scale of the 
shear transformation. 

PACS numbers: 62.20.F-, 61.43.Fs, 83.10.Rs 
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I. INTRODUCTION 


The mechanical properties of bulk metallic glasses, such as high strength and low ductil¬ 
ity, are both of fundamental scientific interest and technological importance jl]. It is now 
well recognized that the plastic deformation of metallic glasses below their glass transition 
temperature involves irreversible rearrangements of small clusters of atoms |2] - The plastic 
flow of amorphous materials in response to applied shear stress can be described using the 
shear transformation zone model, which takes into account the density and internal state of 
the localized zones S. , n recent years, various deformation mechanisms including elemen- 
tary plastic events and shear band formation were studied at different length and time scales 
using atomistic simulations and finite element modeling { 4 ]. Notably, the energy landscape 
analysis had shown that a large strain cycle rejuvenates the glass by increasing the poten¬ 


tial energy, whi 
energy minima 


e a small strain cycle overages the glass by moving the system to deeper 
5]. However, many essential features of the deformation process in strained 
amorphous systems including a correlation between localized plastic events and distribution 
of avalanches are not fully understood. 

The effect of inertia on steadily sheared disordered solids in the athermal quasistatic 
limit was examined in two and three dimensions using molecular dynamics simulations j(3]. 
It was found that the distribution of avalanche sizes obeys a power-law decay over about 
three orders of magnitude in drops of the potential energy density and shear stress, and 
the volume of plastically deformed regions is proportional to the energy dissipated in an 
avalanche j(3]. I 11 the underdamped regime, the system can be carried over successive energy 
barriers to progressively lower minima leading to large avalanches, while in the overdamped 
case, avalanches are smaller and they typically consist of several disconnected regions ori¬ 
ented along diagonal lines js]. It was also shown that at finite strain rates and zero temper¬ 
ature, the correlation between local plastic events remains relevant, and the avalanche size 


scales as the inverse square root of strain rate in two dimensions 


0 - 


In the last few years, a number of studies investigated oscillatory shear response of amor¬ 
phous materials using atomistic simulations js-12], continuum modeling 13], and experi¬ 


mental measurements 


14 


18]. It was found that below a certain strain amplitude, the dis¬ 


ordered systems gradually evolve into dissipative limit cycles and particle rearrangements 
remain reversible, thus retaining memory of the initial conditions |8Hl0i. 12, [l4j]. The num- 
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ber of back and forth cycles required to reach steady state increases as the critical strain 
amplitude is approached from below Js, 10]. Surprisingly, it was shown that the cyclic defor¬ 
mation is accompanied by plastic rearrangements of atoms that are reversed by the end of 


each cycle 


10 


14], With further increasing strain amplitude, particle displacements become 
irreversible leading to a diffusive behavior and structural relaxation BBEQ. 

The elastic response of a two-dimensional amorphous solid to a localized shear trans¬ 
formation was recently studied via molecular dynamics simulations in different damping 


regimes 


19 


20]. In this process, about twenty atoms within a circular inclusion were 


in¬ 


stantaneously sheared in a quiescent system and the time evolution of the displacement held 
was measured. It was demonstrated that the stationary solution for the disorder-averaged 
displacement held has a quadrupolar symmetry and it agrees well with the predictions of 


the continuum elasticity theory [19, 2l|. It was further observed that the transient regime is 


strongly dependent on the damping dynamics and the time dependence of the displacement 
held obtained from molecular dynamics simulations agrees with the continuum solution in 
the overdamped case at large times 19]. The numerical analysis based on the finite element 
method that takes into account microscopic viscosity and the local elastic constants showed 
that the temporal evolution of the disorder-averaged displacement held is similar to the 


propagation of the clastic signal in a uniform medium 
In the previous study 


20 ]. 


22 ], molecular dynamics simulations were carried out to investigate 


the inhuence of a local shear transformation on plastic deformation of a three-dimensional 
model glass. The shear transformation was introduced in a quiescent system via a spherical 
inclusion that was gradually strained into an ellipsoid and then converted back into the 
sphere during a finite time interval. It was demonstrated that at strain amplitudes above a 
few percent, the structural relaxation of the material involved localized plastic events that 
were identihed using the cage detection algorithm 23]. The spatial distribution of clusters 
of cage jumps and their radial density prohles were studied for various damping conditions 
and durations of the shear event. Interestingly, it was found that the density prohles of 


cage jumps are well described by a universal function multiplied 
on the friction coefficient and the shear transformation time scale 


w a factor that depends 


22], It remained unclear, 


however, how this factor varies with the strain amplitude and whether it can be expressed 
as a function of a single variable. 
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In this paper, the plastic response of the amorphous material to a reversible shear trans¬ 


formation is examined over a wide range of damping conditions and oscillation time sea 
The analysis of the density profiles of cage jumps presented in the previous study 22 


es. 

is 


extended further to describe the effect of strain amplitude on the profile shape and the 
dependence of the density maximum on the friction coefficient and oscillation period. In 
particular, it is demonstrated that, at sufficiently slow transformation rates, the peak value 
of the density profiles is a function of the ratio of the friction coefficient and the time scale 
of the shear transformation and that it strongly depends on the strain amplitude. 

The rest of the paper is structured as follows. The details of molecular dynamics simu¬ 
lation model are described in the next section. The analysis of the radial density profiles of 
cage jumps as a function of the time scale of the shear event, friction coefficient, and the 
strain amplitude is presented in Sec. lIIIl The conclusions are provided in the final section. 


II. MOLECULAR DYNAMICS (MD) SIMULATIONS 


The simulated system consists of N = 10 000 particles confined in a three-dimensional cell 


as shown in Fig.lU We used a standard model of a glass-forming Lennard-Jones (LJ) binary 
mixture introduced by Kob and Andersen 24], In this model, any two particles a, j3 = A, B 


interact through the Lennard-Jones (LJ) potential 

V a p{r) = 4e a/3 


12 


U, 


a/3 


( 1 ) 


with the parameters Eaa = 1-0, £ab = 1-5, Ebb = 0.5, (Jab = 0.8, ctbb = 0.88, and 
rriA = rriB 24], The LJ potential is truncated at the cutoff radius r c>a/ j = 2.245 cr,^ 25]. 
The units of length, mass and energy are chosen <r = oaa , m = n^A, and e = Eaa , and, 
consequently, the unit of time is defined r = a^m/E. The simulations were performed at a 
constant density p = Pa + Pb = 1-2 a~ 3 and the linear size of the cubic box is L — 20.27 a. 
Periodic boundary conditions were applied along the x, y, and £ directions. 

The motion of particles is governed by the classical Langevin dynamics. For example, 
the equation of motion in the x direction for the i-th particle of mass m is given by 


mxj 


mTxi = + fi , 




dxi 


( 2 ) 
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where Y is the friction coefficient and /* is a random force with zero mean and vari¬ 
ance ( fi(0)fj{t )) = ^rnksTY5(t)8ij determined by the fluctuation-dissipation theorem 26]. 
The Langevin temperature is fixed T = 10 ~ 2 £/&b, where ks is the Boltzmann constant. 
The equations of motion were integrated numerically using the fifth-order Gear predictor- 
corrector algorithm [27] with a time step A tMD — 0.005 r. Different realizations of disorder 
were prepared by quenching the system from the temperature l.le/ks, which is well above 
T g ~ 0A5e/kB [241], to the final temperature T = 10 ~ 2 e/ks with the rate of 10 ~ 5 e/ksT. 

We next describe the deformation protocol for the reversible shear transformation. The 
inclusion atoms were identified within a sphere of radius vq = 3 a, which is located at the 
center of the simulation cell (see Fig. [I]). The average number of atoms in the inclusion is 
about 135. First, the positions of the inclusion atoms were kept fixed while the system was 
aged for about 500 r at the temperature 10 ~ 2 e/kB- The spherical inclusion was gradually 
strained into an ellipsoid and then converted back into the sphere during the time interval 
r,;. Note that the major axis of the ellipsoid was oriented along one of the diagonals of the 
simulation box, and the volume of the inclusion was kept constant during the transformation. 
In our study, the shear strain is defined as the ratio of the ellipsoid semi-major axis to the 
sphere radius = 3 cr. The variation of stain as a function of time from zero to r* is described 
by the following equation 

e(t) = e 0 sin(7rt/Ti), (3) 


where eo is the strain amplitude and r, is the time scale of the shear event. After the shear 
transformation, the positions of inclusion atoms were kept fixed at their original positions 
and the system was allowed to equilibrate for the additional time interval 10 3 r. This time 
interval is larger than the damping time 1 /Y for the smallest value of the friction coefficient 
T = 0.01 r _1 considered in the present study. Several test simulations were performed 
when the system was equilibrated for a larger time interval of 2 x 10 3 r after the shear 
transformation in order to verify that the results remain unchanged. Finally, the average 
atom positions were computed before and after the shear transformation and then analyzed 
in 768 independent samples. 
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III. RESULTS 


At mechanical equilibrium, the atomic structure of the model glass shows no long-range 
order while each atom remains trapped in a cage composed of its neighbors during the time 
scale of the computer simulation at the studied temperature. The plastic deformation of the 
amorphous material was induced by a reversible shear transformation of a spherical inclusion 
and studied at different damping conditions and time scales of the shear event. In our setup, 
the inclusion atoms were displaced to form an ellipsoid with the major axis parallel to the 
(1,1,1) direction (see Fig.[l]) in order to reduce the effect of periodic boundary conditions. It 
was observed that at sufficiently small strain amplitudes (below a few percent), the 


response is elastic and all atoms return to their cages after the shear transformation 


sy s 


tern 


191 , 


22 ], 


In the present study, the analysis of particle positions was performed in the plastic regime 
when the strain amplitude was varied in the range 0.2 ^ eo ^ 0.4. We fold that at smaller 
values of eo, an accurate analysis of particle displacements requires averaging over larger 
number of independent systems, while at larger strain amplitudes, the relative distance be¬ 
tween inclusion atoms becomes comparable to the molecular diameter, thus creating voids 
at the surface of the inclusion during the shear transformation process. Irreversible rear¬ 
rangements of atoms in the material triggered by the shear transformation were identified 
using the cage detection algorithm 23]- Visual inspection of snapshots of the simulated 
system revealed that cage jumps tend to aggregate into relatively compact clusters, which 
are predominantly located near the inclusion where the deformation of the material during 


the shear transformation is larger 


22 ], 


It was shown in the previous study 


22 ] that the clusters of cage jumps are approximately 


power-law distributed with an exponent that depends on the strain amplitude. In general, 
it is expected that the density of cage jumps will decay away from the center of the inclusion 
because the maximum local strain in the material during the reversible shear transformation 
decreases as a function of the radial distance. For example, the displacement and strain fields 
were calculated analytically for a spherical inclusion that was strained into an ellipsoid in a 


two-dimensional plane, while the third direction reminded neutral 28]. It was shown that 
in a stationary regime, the strain field has a quadrupolar symmetry and it decays as 1/r 3 
from the center of the inclusion in three dimensions 128]. 


Averaged density profiles of cage jumps as a function of the radial distance from the 
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center of the inclusion are plotted in Fig.[2] for the strain amplitude eo = 0.3 in the regime 
of intermediate damping T = 1.0 r _1 . Several important features are evident. First, the 
density of cage jumps is reduced within about two atomic diameters from the surface of 
the inclusion. This effect originates from the reversible motion of the inclusion atoms that 
effectively form a part of a cage for the neighboring atoms of the material, thus reducing 
the probability of their irreversible displacements 22], Second, the density profiles exhibit 
a maximum at r ~ (5 — 6) a and then decay with further increasing radial distance. It was 
previously shown that the rate of decay for r > 6 cr correlates well with the local deformation 
of the material, which was estimated from the relative displacement of neighboring particles 
after the spherical inclusion was irreversibly strained into an ellipsoid 22], 


In general, we find that the average density of cage jumps increases with increasing time 
scale of the shear transformation (e.g., see Fig. [2]). At small values of T; = 5r and 10 r, 
the time scale of the shear event is comparable with the time it takes for sound waves 
to propagate across the system, and thus the probability of formation of large clusters is 
reduced. In contrast, when t* > 50 r, the shear stress from the deforming inclusion can 
induce larger clusters of cage jumps, which in turn might trigger other irreversible events in 
the system. In addition, as the damping rate decreases, the effect of inertia becomes more 
important, leading to larger avalanches during the shear transformation process, and, as a 
result, larger density of cage jumps. This trend was also identified in sheared disordered 
solids in the athermal quasistatic limit by examining the critical scaling of avalanches at 
different damping rates [6]. 


It was further noticed in the previous study 22] that the density profiles of cage jumps 
for different values of r* and T can be made to collapse onto a master curve if p{r ) is 
divided by a scaling factor. In Fig.[3l we plot the average density profiles normalized by 
the corresponding density peak p m for different values of V and strain amplitudes e 0 = 0.2, 
0.3 and 0.4. It is evident that for each value of the strain amplitude, the rescaled density 
profiles p(r)/p m collapse on the master curves. Note that the data in Fig. [3] for the two 
lower values eo = 0.2 and 0.3 are shifted for clarity. It can be seen that the location of the 
maximum of p{r)/p m and the slope of decay for r > 6 a depend of the strain amplitude. 
The deviation from the 1/r 3 dependence, which describes the decay of the local strain away 


from an elliptical inclusion 


, for larger strain amplitudes eo = 0.3 and 0.4 might be due 
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to the finite system size. Remember that at r>L/2 ^ 10 cr some atoms interact with their 
neighbors via periodic boundary conditions but the local strain is in general not the same 
across periodic boundaries during the shear transformation. 

In our study, the maximum of the density profiles of cage jumps p m was estimated in 
a wide range of parameter values, i.e., 0.01 ^ Tr ^ 10 and 5 r ^ T; 10 3 r. Figured 
shows a contour plot of the density peaks p m as a function of the friction coefficient and 
the shear transformation time scale for the strain amplitude eo = 0.3. It can be seen that 
the density landscape is quite complex, but the trends are clear. Namely, at small r,: in the 
overdamped regime, the deformation of material during the shear transformation is minimal, 
and thus the density of cage jumps is relatively small. In the opposite limit, when the shear 
transformation is very slow and the system dynamics is underdamped, the deformation of 
material is largest, facilitating the formation of large clusters, and the density of cage jumps 
saturates to a maximum value of about 0.037 c~ 3 (see Fig. [4]). 

It can be further observed that the contour lines in Fig. 0] approximately follow a linear 
dependence between the friction coefficient and the time scale of the shear event (see the 
straight line with unit slope in Fig.0|). This correlation holds over the whole range of 
parameters Y and r*, except for the data points within the dashed region in Fig. 01 This, 
in turn, implies that the ratio T/t; computed along a contour level will correspond to the 
same value of the density peak. Therefore, it is expected that the data reported outside of 
the dashed region in Fig. 0] can be collapsed onto a master curve if replotted as a function 
of T/rj. Another argument for using T/t; is that the time dependence of the continuum 
displacement field after an instantaneous shear transformation in the overdamped regime is 
roughly proportional to the factor e~ r2r ^ flt , where p is the shear modulus [jij. Therefore, it 
follows that the local displacement field, which can trigger an irreversible rearrangement of 
atoms at a distance r from an inclusion, depends on the ratio Y/t. 

The density peaks of cage jumps are first plotted in Fig. [5] as a function of the ratio Y/t % 
for the strain amplitudes eo = 0.2, 0.3 and 0.4 and the same range of parameter values as in 
Fig. 01 The collapse of the data on three distinct master curves is satisfactory; however, the 
data are somewhat scattered at intermediate values of Y/t % . As anticipated, the scattered 
data in Fig. [5] were evaluated for the parameters Y and r t within the dashed region in Fig. 01 
Without these data points, the dependence of the density peaks on the ratio T/t* is shown in 


Fig. [6] It is evident that for T /r \ < 0.01, the density peak saturates to a constant value that 
corresponds to the largest plastic deformation for each strain amplitude. In contrast, with 
increasing T/t*, the density peak gradually crosses over to a power-law decay as a function 
of T/Tj. In this regime, the density of cage jumps is reduced due to either small time scale 
Tj or large friction coefficient. 

It is apparent that the shape of the curves shown in Fig.[6]is very similar, suggesting that 
they might be different by a factor that depends on the strain amplitude. Indeed, we found 
that when the density peaks for each strain amplitude are divided by e ( f, the data collapse 
onto a single master curve (see Fig. ED- The resulting master curve extends over about five 
orders of magnitude in T/t*. Notice that the slope of the decay at T/t* > 0.01 depends 
slightly on the strain amplitude. The existence of the plateau in Fig. [7] implies that the 
largest value of p m can be obtained in the limiting case of a very slow shear transformation 
for any damping conditions, i.e., when T/t* —> 0. The significance of the value 5 for the 
exponent is at present not clear. We note, however, that the exponent was estimated based 
only on three data points for eo and the critical value of the strain amplitude that marks 
the onset of irreversible deformation was not determined in our study. The relatively strong 
dependence of p m on the strain amplitude might be due the quadrupolar symmetry of the 
strain held that can trigger increasingly large clusters of cage jumps around the inclusion 
upon increasing strain amplitude. 


IV. CONCLUSIONS 

In this paper, we have examined the structural relaxation in a three-dimensional amor¬ 
phous material induced by a reversible shear transformation using molecular dynamics simu¬ 
lations. The material was deformed by straining a spherical inclusion into an ellipsoid of the 
same volume and then converting it back into the sphere. We found that at sufficiently large 
strain amplitude of the shear transformation, some particles undergo irreversible displace¬ 
ments that were identified using the cage detection algorithm. The density profiles of cage 
jumps exhibit a distinct maximum near the surface of the inclusion followed by a power-law 
decay as a function of the radial distance. At a given strain amplitude, the density profiles 
are self-similar when scaled by the density maximum, which in turn depends on the damping 
rate and duration of the shear transformation. Moreover, it was demonstrated that the data 
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for the peak value of the density profiles can be collapsed onto a master curve when plotted 
as a function of the ratio of the friction coefficient and the oscillation period. Overall, these 
findings indicate that the density of cage jumps around the inclusion becomes larger in the 
cases of weakly damped dynamics or slow shear transformation. 
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FIG. 1: (Color online) A snapshot of the instantaneous configuration of atoms of type A (large blue 
circles) and type B (small red circles) in the binary (80:20) LJ mixture. The spherical inclusion 
is located at the center of the periodic cell (black circle). The reversible shear transformation is 
applied to the inclusion atoms, which are gradually strained into an ellipsoid of the same volume 
(dashed ellipse) and then returned to their original positions. 
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FIG. 2: (Color online) Averaged radial density profiles of cage jumps for the friction coefficient 
r = 1.0 r” 1 and the strain amplitude eo = 0.3. The time scale of the shear transformation is 
Tj/r = 5, 10, 50, 100, 300, 500, and 1000 from bottom to top. 
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FIG. 3: (Color online) Log-log plot of the rescaled density profiles of cage jumps for the strain 
amplitudes eo = 0.2, 0.3, and 0.4 and the shear transformation time scale r,- = 100 r. For each 
strain amplitude, the friction coefficient is Tr = 0.01 (solid black curve), 0.1 (dashed blue curve), 
1.0 (dash-dotted red curve), and 10 (solid green curve). The data for eo = 0.2 and 0.3 are displaced 
vertically for clarity. The straight solid and dashed lines indicate slopes —3 and —4 respectively. 
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FIG. 4: (Color online) A contour plot of the density peaks p m (in units of <r~ 3 ) as a function of the 
friction coefficient T and the time scale of the shear event Tj for the strain amplitude eo = 0.3. The 
contour levels are specified in the legend. Open square symbols indicate individual data points. 
The straight line with unit slope is shown as a reference. The data points within the dashed region 
were excluded in the analysis of p rn presented in Figs. [6] and [7] (see text for details). 
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FIG. 5: (Color online) The variation of the density peak p m as a function of the ratio T/t* for the 
strain amplitudes eo = 0.2 (o), 0.3 (o), and 0.4 (A). The friction coefficient and the time scale of 
the shear event vary in the ranges 0.01 ^ Tr ^ 10 and 5r ^ r* ^ 10 3 r, respectively. Error bars 
are about the size of the symbols. 
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FIG. 6: (Color online) The variation of the density peak p m as a function of the ratio T/r* for the 
strain amplitudes eo = 0.2 (o), 0.3 (o), and 0.4 (A). The same data as in Fig.[5] except for Tr ^ 1 
and Tj ^ 10 t (i.e., except for the data points within the dashed region shown in Fig.0|). 
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FIG. 7: (Color online) Master plot of p m /e q versus T/t,. The same data as in Fig.[6l The straight 
line with a slope —0.5 is shown for reference. 
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